from math import cos, floor, pi

from numpy.core import shape_base
from numpy.core.fromnumeric import size
from numpy.lib.npyio import save
from gear import Gear, GearInit, PlanetaryGear, PlanetInit
import numpy as np
from matplotlib import pyplot as plt
import shichuan_stiff as st
import os

from matplotlib import rcParams
#rcParams.update({'font.size': 11, 'text.usetex': True})
rcParams.update({'font.size': 13, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'})

def H(e):
    if e >= 0:
        return 1
    else:
        return 0


def Ul_h_i(e, i, epsilon_a, epsilon_b):

    if epsilon_b > 0:
        return (1/epsilon_b) * (
            H(e-i)*(e-i) - H(e-epsilon_b-i)*(e-epsilon_b-i) - H(e-epsilon_a-i)*(e-epsilon_a-i) +
            H(e-(epsilon_a+epsilon_b)-i)*(e-(epsilon_a+epsilon_b)-i)
        )
    else:
        return H(e-i)-H(e-epsilon_a-i)


def Tl_h(e, epsilon_a, epsilon_b):
    return H(e) - H(e-(epsilon_a+epsilon_b))


def L_h(e, epsilon_a, epsilon_b):
    tmp = 0
    i = -floor(epsilon_a+epsilon_b)
    while (i <= floor(epsilon_a+epsilon_b)):
        tmp += Ul_h_i(e, i, epsilon_a, epsilon_b)
        i += 1
    return tmp * Tl_h(e, epsilon_a, epsilon_b)


def Kmax(gg: Gear, methods='iso'):
    if methods == 'iso':
        zv1 = gg.z1/cos(gg.beta)**3  # 当量齿数
        zv2 = gg.z2/cos(gg.beta)**3
        x1 = gg.x1
        x2 = gg.x2
        alpha = gg.alpha
        betab = gg.betab

        q = (
            0.04723+0.15551/zv1+0.25791/zv2-0.00635*x1-0.11654*x1/zv1 -
            0.00193*x2-0.24188*x2/zv2+0.00529*x1**2+0.00182*x2**2
        )

        return 1/q * cos(alpha) * cos(gg.beta)
    else:
        tmp, _ = st.stiffness(gg)
        # print(tmp[0][2])
        return max(list(tmp[0][2]))


def k_u_i(e, i, epsilon_a, epsilon_b, ak=0.8):
    return 4*(ak-1)*(e-i)**2/(epsilon_a+epsilon_b)**2 - 4*(ak-1)*(e-i)/(epsilon_a+epsilon_b) + ak


def K_u_i_Devide_Kmax(e, i, epsilon_a, epsilon_b, b, betab, ak=0.8):
    return k_u_i(e, i, epsilon_a, epsilon_b, ak)*cos(betab)/b


def KL_h(e, epsilon_a, epsilon_b, b, betab, ak=0.8):
    tmp = 0
    i = -floor(epsilon_a+epsilon_b)
    while (i <= floor(epsilon_a+epsilon_b)):
        tmp += Ul_h_i(e, i, epsilon_a, epsilon_b)*b/cos(betab) * \
            K_u_i_Devide_Kmax(e, i, epsilon_a, epsilon_b, b, betab, ak)
        i += 1
    return tmp * Tl_h(e, epsilon_a, epsilon_b)


def stiffness_hs(gg: Gear, xstep=0.01, savefilename=None, show=False, kmaxcalc=True, ak=0.8, kmaxmethod='iso'):
    '''
    return 斜齿轮or直齿轮的时变啮合刚度
    '''
    # gg.calc()
    epsilon_a = gg.ContactRatiot
    epsilon_b = gg.ContactRatioBeta
    print(epsilon_a, epsilon_b, epsilon_a/epsilon_b)
    b = min([gg.b1, gg.b2])
    betab = gg.betab
    elist = np.arange(0, gg.ContactRatio, xstep)

    if kmaxcalc:
        kmax = Kmax(gg, kmaxmethod)
        # print(kmax)
    else:
        kmax = 1  # 计算无量纲刚度
    KallList = []
    KsingleList_0 = []
    KsingleList_p1 = []
    KsingleList_1 = []
    for e in elist:
        KallList.append(KL_h(e, epsilon_a, epsilon_b, b, betab, ak=ak)*kmax)
        KsingleList_0.append(k_u_i(e, 0, epsilon_a, epsilon_b, ak=ak)*kmax *
                             Ul_h_i(e, 0, epsilon_a, epsilon_b)*Tl_h(e, epsilon_a, epsilon_b))
        KsingleList_p1.append(k_u_i(e, -1, epsilon_a, epsilon_b, ak=ak)*kmax *
                              Ul_h_i(e, -1, epsilon_a, epsilon_b)*Tl_h(e, epsilon_a, epsilon_b))
        KsingleList_1.append(k_u_i(e, 1, epsilon_a, epsilon_b, ak=ak)*kmax *
                             Ul_h_i(e, 1, epsilon_a, epsilon_b)*Tl_h(e, epsilon_a, epsilon_b))

    if (savefilename != None) or show:  # 如果既不储存也不显示，就不作图，仅仅返数据
        plt.figure(figsize=(7, 7))
        plt.plot(elist, KallList, 'r-', label='K_all')
        plt.plot(elist, KsingleList_0, 'b--', label='K(i=0)')
        plt.plot(elist, KsingleList_p1, 'g--', label='K(i=-1)')
        plt.plot(elist, KsingleList_1, 'g--', label='K(i=1)')
        plt.xlabel('${\\xi}$')
        if kmaxcalc:
            plt.ylabel('${K^{LS}}$/(N/(mm${\\cdot\\mu}$m))')
            plt.title('Stiffness(QSEMmodel):${K^{LS}}$ -- ${\\xi}$')
        else:
            plt.ylabel('${K^{LS}/K_{\max}}$/(N/(mm${\\cdot\\mu}$m))')
            plt.title('Stiffness(QSEMmodel):${K^{LS}/K_{\max}}$ -- ${\\xi}$')
        plt.grid()
        plt.legend()
        if savefilename != None:
            if not os.path.exists('./pic'):
                os.mkdir('./pic')
            plt.savefig("./pic/"+savefilename, dpi=400, bbox_inches='tight')
        if show:
            plt.show()
    pass

    return elist, KallList, KsingleList_0, KsingleList_p1, KsingleList_1


def comparespur(gg: Gear):

    e1, k1, ks0, ksp1, ks1 = stiffness_hs(
        gg, xstep=0.005, kmaxmethod='shichuan')
    xs2, _ = st.stiffness(gg, xstep=0.005, savefilename=None, show=False)

    plt.figure(figsize=(14, 7))
    plt.subplot(1, 2, 1)
    plt.plot(e1, k1, 'r-', label='sum')
    plt.plot(e1, ks0, 'b--', label='i=0')
    plt.plot(e1, ksp1, 'g--', label='i=-1')
    plt.plot(e1, ks1, 'g--', label='i=1')
    plt.xlabel('${\\xi}$')
    plt.ylabel('K/(N/(mm${\\cdot\\mu}$m))')
    plt.title('K-${\\xi}$')
    plt.ylim([0, 25])
    plt.grid()
    plt.legend()

    plt.subplot(1, 2, 2)
    plt.plot(xs2[0][0]*180/pi, xs2[0][1], 'r-',
             label='sum')
    plt.plot(xs2[0][0]*180/pi, xs2[0][2], 'b--',
             label='i=0')
    #x_range_1and3 = np.append(xrange1, xrange3)
    plt.plot(xs2[1][0]*180/pi, xs2[1][1], 'g--',
             label='i=-1')

    plt.plot(xs2[2][0]*180/pi, xs2[2][1], 'g--',
             label='i=1')
    plt.xlabel('Angle Rotate/(deg)')
    plt.ylabel('K/(N/(mm${\\cdot\\mu}$m))')
    plt.title('K-Angle')
    plt.ylim([0, 22])
    plt.grid()
    plt.legend()
    plt.savefig("./pic/"+'comparespur.png', dpi=400, bbox_inches='tight')
    plt.show()


def changebeta():
    initg = GearInit(b=(40, 40),)
    plt.figure(figsize=(7, 5))
    typelist = ['-', '--', ':', '-.']
    i = 0
    for beta_deg in np.arange(14, 21, 2):
        gg = Gear(4.5, 16, 24, 0.1817, 0.1715, initg, beta_deg=beta_deg)
        gg.calc()
        e1, k1, ks0, ksp1, ks1 = stiffness_hs(
            gg, xstep=0.005, kmaxmethod='iso', kmaxcalc=False)
        char = '${\\beta}$=%d' % (beta_deg)
        plt.plot(e1, k1, typelist[i], label=char)
        i += 1
    plt.xlabel('${\\xi}$')
    plt.ylabel('K/(N/(mm${\\cdot\\mu}$m))')
    plt.title('K-${\\xi}$')
    plt.legend()
    plt.grid()
    # plt.ylim([10,25])
    #plt.savefig("./pic/"+'QSMEgear1change.png', dpi=400, bbox_inches='tight')
    plt.show()


def testepdiv():

    typelist = ['-', '--', ':', '-.']
    colorlist = ['red', 'orange', 'green', 'blue', 'black']
    ilist = np.arange(0.8, 1.2+0.01, 0.1)
    Kmat = []
    plt.figure(figsize=(7, 5))
    ak = 0.8
    j = int(len(ilist)/2)
    for i in ilist:
        epa = 1.4
        epb = epa * (i/epa)
        elist = np.arange(0, epa+epb, 0.005)
        KallList = []
        KsingleList_0 = []
        KsingleList_p1 = []
        KsingleList_1 = []
        for e in elist:
            KallList.append(KL_h(e, epa, epb, 2, 0, ak=ak))
            '''
            KsingleList_0.append(k_u_i(e, 0, epa, epb, ak=ak)*
                                Ul_h_i(e, 0, epa, epb)*Tl_h(e, epa, epb))
            KsingleList_p1.append(k_u_i(e, -1, epa, epb, ak=ak) *
                                Ul_h_i(e, -1, epa, epb)*Tl_h(e, epa, epb))
            KsingleList_1.append(k_u_i(e, 1, epa, epb, ak=ak)*
                                Ul_h_i(e, 1, epa, epb)*Tl_h(e, epa, epb))
            '''
        Kmat.append(KallList)

        plt.plot(elist, KallList, typelist[j % len(typelist)], color=colorlist[j % len(
            colorlist)], label='$\\frac{ep_{a}}{ep_{b}}=%.1fep_{a}$' % (1/epb))

        plt.annotate('$%.1f$' % (
            1/epb), xy=(elist[-1]-0.05, KallList[-1]-0.03), color=colorlist[j % len(colorlist)], size=9)
        j += 1

    '''
    plt.plot(elist, KsingleList_0, 'b--', label='K(i=0)')
    plt.plot(elist, KsingleList_p1, 'g--', label='K(i=-1)')
    plt.plot(elist, KsingleList_1, 'g--', label='K(i=1)')
    '''

    plt.xlabel('${\\xi}$')
    plt.ylabel('${K^{LS}/K_{\max}}$/(N/(mm${\\cdot\\mu}$m))')
    plt.title(
        'Stiffness(QSEMmodel):${K^{LS}/K_{\max}}$ -- ${\\xi};ep_{a}=%.2f$' % (epa))
    plt.ylim([epa-0.5, epa+0.2])
    plt.grid()
    plt.legend()

    #plt.savefig("./pic/"+'epdiv%d.png'%(epa*10), dpi=400, bbox_inches='tight')
    plt.show()


def testFixedall():
    '''
    测试总重合度固定
    '''
    eall = 3.5

    epalist = np.arange(1.1, 1.9+0.01, 0.2)

    typelist = ['-', '--', ':', '-.']
    colorlist = ['red', 'orange', 'green', 'blue', 'black']
    Kmat = []
    plt.figure(figsize=(7, 5))
    ak = 0.8
    j=0

    for epa in epalist:
        epb = eall - epa
        elist = np.arange(0, epa+epb, 0.005)
        KallList = []
        for e in elist:
            KallList.append(KL_h(e, epa, epb, 2, 0, ak=ak))

        Kmat.append(KallList)

        plt.plot(elist, KallList, typelist[j % len(typelist)], color=colorlist[j % len(
            colorlist)], label='$ep_{a}=%.1f;ep_{b}=%.1f$' % (epa,epb))

        plt.annotate('$ep_{b}=%.1f$' % (
            epb), xy=(elist[-1]-0.05, KallList[-1]-0.03), color=colorlist[j % len(colorlist)], size=9)
        j += 1

    '''
    plt.plot(elist, KsingleList_0, 'b--', label='K(i=0)')
    plt.plot(elist, KsingleList_p1, 'g--', label='K(i=-1)')
    plt.plot(elist, KsingleList_1, 'g--', label='K(i=1)')
    '''

    plt.xlabel('${\\xi}$')
    plt.ylabel('${K^{LS}/K_{\max}}$/(N/(mm${\\cdot\\mu}$m))')
    plt.title(
        'Stiffness(QSEMmodel):${K^{LS}/K_{\max}}$ -- ${\\xi};ep_{all}=%.2f$' % (eall))
    #plt.ylim([epa-0.5, epa+0.2])
    plt.grid()
    plt.legend()

    plt.savefig("./pic/"+'epFixed_Epall%d.png'%(eall*10), dpi=400, bbox_inches='tight')
    plt.show()

def articledraw():
    #plt.rc('font',family='Times New Roman')
    
    epa = 1.4
    epb = 0.6
    '''
    elist = np.arange(-1.6, 3, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,0))
        ss0.append(Ul_h_i(e,0,epa,0))
        ss1.append(Ul_h_i(e,1,epa,0))
    
    plt.figure(figsize=(4,4))
    plt.plot(elist,ssp1,'--',color='black',label='$i=-1$')
    plt.plot(elist,ss0,'-',color='black',label='$i=0$')
    plt.plot(elist,ss1,':',color='black',label='$i=1$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$s^{s}(\\xi)$')
    plt.legend(loc='lower right')
    plt.grid()
    plt.savefig("./picar/"+'ss.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    elist = np.arange(-1, 3, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb))
    
    plt.figure(figsize=(4,4))
    plt.plot(elist,ssp1,'--',color='black',label='$i=-1$')
    plt.plot(elist,ss0,'-',color='black',label='$i=0$')
    plt.plot(elist,ss1,':',color='black',label='$i=1$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$s^{h}(\\xi)$')
    plt.legend(loc='lower right')
    plt.grid()
    plt.savefig("./picar/"+'sh.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    elist = np.arange(-1, 3, 0.005)
    ss0 = []
    for e in elist:
        ss0.append(Tl_h(e,epa,epb))
    
    plt.figure(figsize=(4,4))
    plt.plot(elist,ss0,'-',color='black',label='$i=0$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$T(\\xi)$')
    plt.legend(loc='upper right')
    plt.grid()
    plt.savefig("./picar/"+'T.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    elist = np.arange(-0.5, 2.5, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    ssall = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,epb)*Tl_h(e,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb)*Tl_h(e,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb)*Tl_h(e,epa,epb))
        ssall.append(L_h(e,epa,epb))
    
    plt.figure(figsize=(4,4))
    plt.plot(elist,ssp1,'--',color='black',label='$i=-1$')
    plt.plot(elist,ss0,'-.',color='black',label='$i=0$')
    plt.plot(elist,ss1,':',color='black',label='$i=1$')
    plt.plot(elist,ssall,'-',color='black',label='$S(\\xi)$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$S(\\xi)$')
    plt.legend(loc='lower right')
    plt.grid()
    plt.savefig("./picar/"+'Sall.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    elist = np.arange(0, 2, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    ssall = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,epb)*Tl_h(e,epa,epb)*k_u_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb)*Tl_h(e,epa,epb)*k_u_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb)*Tl_h(e,epa,epb)*k_u_i(e,1,epa,epb))
        ssall.append(KL_h(e,epa,epb,1,0))
    
    plt.figure(figsize=(4,4))
    plt.plot(elist,ssp1,'--',color='black',label='$k_{-1}(\\xi)$')
    plt.plot(elist,ss0,'-.',color='black',label='$k_{0}(\\xi)$')
    plt.plot(elist,ss1,':',color='black',label='$k_{1}(\\xi)$')
    plt.plot(elist,ssall,'-',color='black',label='$K(\\xi)$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='lower right')
    plt.grid()
    plt.savefig("./picar/"+'Kall.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    ##2
    epa=1.4
    epb=0.6
    elist = np.arange(0, epa+epb, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb))
    
    plt.figure(figsize=(4,4))
    #plt.plot(elist,ssp1,'--',color='black',label='i=-1')
    plt.plot(elist,ss0,'-',color='black',label='$i=0$')
    plt.plot([epb,epb],[-0.05,1.05],'--',color='black')
    plt.plot([epa,epa],[-0.05,1.05],'--',color='black')
    plt.text(epb+0.02,-0.05,'$\\varepsilon_{b}$')
    plt.text(epa+0.02,-0.05,'$\\varepsilon_{a}$')
    #plt.plot(elist,ss1,':',color='black',label='i=1')
    plt.xlabel('$\\xi$')
    plt.ylabel('$s(\\xi)$')
    plt.legend(loc='upper right')
    plt.grid()
    plt.savefig("./picar/"+'epbepa.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    epa=0.8
    epb=1.3
    elist = np.arange(0, epa+epb, 0.005)
    ssp1 = []
    ss0 = []
    ss1 = []
    for e in elist:
        ssp1.append(Ul_h_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb))
    
    plt.figure(figsize=(4,4))
    #plt.plot(elist,ssp1,'--',color='black',label='i=-1')
    plt.plot(elist,ss0,'-',color='black',label='$i=0$')
    plt.plot([epb,epb],[-0.05,1.05],'--',color='black')
    plt.plot([epa,epa],[-0.05,1.05],'--',color='black')
    plt.text(epb+0.02,-0.05,'$\\varepsilon_{b}$')
    plt.text(epa+0.02,-0.05,'$\\varepsilon_{a}$')
    #plt.plot(elist,ss1,':',color='black',label='i=1')
    plt.xlabel('$\\xi$')
    plt.ylabel('$s(\\xi)$')
    plt.legend(loc='upper right')
    plt.grid()
    plt.savefig("./picar/"+'epaepb.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    epa=1.5
    epb=1
    
    elist = np.arange(-2, 4.5, 0.005)
    ssp2 = []
    ssp1 = []
    ss0 = []
    ss1 = []
    ss2 = []
    ssall = []
    for e in elist:
        ssp2.append(Ul_h_i(e,-2,epa,epb))
        ssp1.append(Ul_h_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb))
        ss2.append(Ul_h_i(e,2,epa,epb))
        ssall.append(L_h(e,epa,epb))
    
    plt.figure(figsize=(5,4))
    xlist=np.arange(-2,2.1,1)
    for x in xlist:
        plt.plot([x,x],[0,-0.05],color='black')
        if x!=2:
            plt.text(x+0.3,-0.05,'$inc↑$')
    
    
    xlist=np.arange(-0.5,3.51,1)
    for x in xlist:
        plt.plot([x,x],[1,1.05],color='black')
        if x!=3.5:
            plt.text(x+0.25,1.02,'$dec↓$')
    
    plt.plot(elist,ssp2,'-',color='gray',label='$s_{-2}(\\xi)$')
    plt.plot(elist,ssp1,'--',color='black',label='$s_{-1}(\\xi)$')
    plt.plot(elist,ss0,'-',color='black',label='$s_{0}(\\xi)$')
    plt.plot(elist,ss1,':',color='black',label='$s_{1}(\\xi)$')
    plt.plot(elist,ss2,'-.',color='black',label='$s_{2}(\\xi)$')
    #plt.plot(elist,ssall,'-',color='black',label='$K(\\xi)$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$S(\\xi)$')
    plt.legend(loc='lower right',fontsize=10)
    plt.grid()
    plt.savefig("./picar/"+'ebnpt.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    epa=1.5
    epb=1
    
    elist = np.arange(0, 2.5, 0.005)
    ssp2 = []
    ssp1 = []
    ss0 = []
    ss1 = []
    ss2 = []
    ssall = []
    for e in elist:
        ssp2.append(Ul_h_i(e,-2,epa,epb))
        ssp1.append(Ul_h_i(e,-1,epa,epb))
        ss0.append(Ul_h_i(e,0,epa,epb))
        ss1.append(Ul_h_i(e,1,epa,epb))
        ss2.append(Ul_h_i(e,2,epa,epb))
        ssall.append(L_h(e,epa,epb))
    
    plt.figure(figsize=(4,4))
    
    # xlist=np.arange(-2,2.1,1)
    # for x in xlist:
    #     plt.plot([x,x],[0,-0.05],color='black')
    #     if x!=2:
    #         plt.text(x+0.3,-0.05,'$inc↑$')
    
    
    # xlist=np.arange(-0.5,3.51,1)
    # for x in xlist:
    #     plt.plot([x,x],[1,1.05],color='black')
    #     if x!=3.5:
    #         plt.text(x+0.25,1.02,'$dec↓$')
    
    #plt.plot(elist,ssp2,'-',color='gray',label='$k_{-2}(\\xi)$')
    plt.plot(elist,ssp1,'--',color='black',label='$s_{-1}(\\xi)$')
    plt.plot(elist,ss0,'-.',color='black',label='$s_{0}(\\xi)$')
    plt.plot(elist,ss1,':',color='black',label='$s_{1}(\\xi)$')
    #plt.plot(elist,ss2,'-.',color='black',label='$k_{2}(\\xi)$')
    plt.plot(elist,ssall,'-',color='black',label='$S(\\xi)$')
    plt.xlabel('$\\xi$')
    plt.ylabel('$S(\\xi)$')
    plt.legend(loc='lower right',fontsize=10)
    plt.grid()
    plt.savefig("./picar/"+'ebnptall.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    epa=0.8
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epb in [0.8,0.9,1,1.1,1.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{b}=%.2f$'%(epb))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='upper right',fontsize=10)
    plt.ylim([0.5,1.2])
    plt.grid()
    plt.savefig("./picar/"+'epa08.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    epa=1.4
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epb in [0.8,0.9,1,1.1,1.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{b}=%.2f$'%(epb))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='upper right',fontsize=10)
    plt.ylim([1.0,1.8])
    plt.grid()
    plt.savefig("./picar/"+'epa14.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    epa=1.7
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epb in [0.8,0.9,1,1.1,1.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{b}=%.2f$'%(epb))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='upper right',fontsize=10)
    plt.ylim([1.3,2.1])
    plt.grid()
    plt.savefig("./picar/"+'epa17.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    epa=1.5
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epb in [0.8,0.9,1,1.1,1.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{b}=%.2f$'%(epb))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='upper right',fontsize=10)
    plt.ylim([1.1,2.0])
    plt.grid()
    plt.savefig("./picar/"+'epa15epb1.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    
    epa=1.5
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epb in [1.8,1.9,2,2.1,2.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{b}=%.2f$'%(epb))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='upper right',fontsize=10)
    plt.ylim([1.25,1.75])
    plt.grid()
    plt.savefig("./picar/"+'epa15epb2.png', dpi=400, bbox_inches='tight')
    plt.show()
    '''
    epb=1.3
    plt.figure(figsize=(5,4))
    typelist=['-','--','-',':','-.']
    colorlist=['gray','black','black','black','black']
    j=0
    for epa in [0.8,0.9,1,1.1,1.2]:
        
        ss1 = []
        elist = np.arange(0, epa+epb, 0.005)
        
        for e in elist:
            ss1.append(KL_h(e,epa,epb,1,0))
        # if j==0:
        #     j+=1
        #     continue
        plt.plot(elist,ss1,typelist[j],color=colorlist[j],label='$\\varepsilon_{a}=%.2f$'%(epa))
        j+=1
        
    #plt.title('$\\varepsilon_{a}=%.2f$'%epa)
    plt.xlabel('$\\xi$')
    plt.ylabel('$K(\\xi)$')
    plt.legend(loc='lower right',fontsize=10)
    plt.ylim([0.5,1.3])
    plt.grid()
    plt.savefig("./picar/"+'epb13epa1.png', dpi=400, bbox_inches='tight')
    plt.show()
    
    pass

if __name__ == '__main__':
    #testepdiv()
    #testFixedall()
    articledraw()
    
    
    
    
    
    '''
    init = PlanetInit()
    pt = PlanetaryGear(1.75, 18, 45, 108, 0.393,
                       0.0285, 0.45, init, beta_deg=20)
    g13 = pt.gear_ac

    initg = GearInit(b=(40, 40),)
    initg.E1 = g13.InitStruct.E1
    initg.E2 = g13.InitStruct.E2
    initg.v1 = g13.InitStruct.v1
    initg.v2 = g13.InitStruct.v2
    
    gg = Gear(4.5, 16, 24, 0.1817, 0.1715, initg, beta_deg=15)
    gg.calc()
    '''
    # comparespur(gg)
    # stiffness_hs(gg,show=True,savefilename='QSMEstiff1.png')
    #init2 = GearInit(b=(40, 40), types='nei')
    #gg2 = Gear(1.75, 45, 108, 0.0285, 0.45, init2, beta_deg=0)
    # gg2.calc()
    # changebeta()
    # stiffness_hs(gg,show=True,kmaxmethod='iso',savefilename='QSMEgear1.png',kmaxcalc=True)
    # st.stiffness(gg2,show=True)

    '''
    epa = 1.6
    epb = 0.9
    plt.figure(figsize=(7,6))
    elist=np.arange(0,epa+epb,0.005)
    ulist=[]
    mlistp1=[]
    mlist0=[]
    mlist1=[]
    ccp1=[]
    cc=[]
    cc1=[]
    klistp1=[]
    klist=[]
    klist1=[]
    ak=0.8
    for e in elist:
        
        klistp1.append(k_u_i(e,-1,epa,epb,0.8)*Tl_h(e, epa, epb))
        klist.append(k_u_i(e,0,epa,epb,0.8)*Tl_h(e, epa, epb))
        klist1.append(k_u_i(e,1,epa,epb,0.8)*Tl_h(e, epa, epb))
        
        mlistp1.append(k_u_i(e, -1, epa, epb, ak=ak)*Ul_h_i(e, -1, epa, epb)*Tl_h(e, epa, epb))
        mlist0.append(k_u_i(e, 0, epa, epb, ak=ak)*Ul_h_i(e, 0, epa, epb)*Tl_h(e, epa, epb))
        mlist1.append(k_u_i(e, 1, epa, epb, ak=ak)*Ul_h_i(e, 1, epa, epb)*Tl_h(e, epa, epb))
        
        ccp1.append(Ul_h_i(e, -1, epa, epb))#*Tl_h(e, epa, epb))
        cc.append(Ul_h_i(e, 0, epa, epb))#*Tl_h(e, epa, epb))
        cc1.append(Ul_h_i(e, 1, epa, epb))#*Tl_h(e, epa, epb))
        
        ulist.append(KL_h(e,epa,epb,2,0,ak))
    
    # plt.plot(elist,klistp1,'r-',label='gear pre')
    # plt.plot(elist,klist,'g-',label='gear now')
    # plt.plot(elist,klist1,'b-',label='gear next')
    # #plt.annotate('Single', xy=(1,1.03), color='blue', size=10)
    
    # plt.plot(elist,ccp1,'r:')
    # plt.plot(elist,cc,'g:')
    # plt.plot(elist,cc1,'b:')
    
    #plt.annotate('ContactLine', xy=(0.5,1.03), color='red', size=10)
    
    # plt.plot(elist,mlistp1,'r--',label='gear pre')
    # plt.plot(elist,mlist0,'g--',label='gear now')
    # # #plt.annotate('Exact Stiffness', xy=(0.7,0.95), color='black', size=10)
    # plt.plot(elist,mlist1,'b--',label='gear next')
    
    color = 'red'
    plt.plot(elist,ulist,'-',color=color)
    plt.annotate('%.2f'%(epb), xy=(elist[-1]+0.01,ulist[-1]), color=color, size=10)
    
    color = 'orange'
    plt.plot(elist,ulist,'-',color=color)
    plt.annotate('%.2f'%(epb), xy=(elist[-1]+0.01,ulist[-1]), color=color, size=10)
    
    
    plt.xlabel('${\\xi}$')
    plt.ylabel('${N(\\xi)}$')
    #plt.title('Contact Line Length: $\\frac{ep_{a}}{ep_{b}}=%.1fep_{a}$'%(1/epb))
    plt.grid()
    #plt.xlim([-1,3])
    plt.ylim([1.1,1.8])
    #plt.xlim([0,2])
    #plt.legend(loc='upper right')
    plt.savefig("./pic/"+'ep2.png', dpi=400, bbox_inches='tight')
    plt.show()
    '''
    
    
    
    '''
    typelist = ['-', '--', ':', '-.']
    typelist = ['-']
    colorlist = ['red', 'orange', 'green', 'blue', 'black']
    ilist = np.arange(1.2+1,1+0.8-0.01, -0.1)
    Kmat = []
    plt.figure(figsize=(7, 5))
    epa = 1.4
    plt.ylim([epa-0.3, epa+0.1])
    #plt.grid()
    ak = 0.8
    j = int(len(ilist)/2)
    for i in ilist:
        
        epb = epa * (i/epa)
        elist = np.arange(0, epa+epb, 0.005)
        KallList = []
        for e in elist:
            KallList.append(KL_h(e, epa, epb, 2, 0, ak=ak))
            
        Kmat.append(KallList)

        plt.plot(elist, KallList, typelist[j % len(typelist)], color=colorlist[j % len(
            colorlist)], label='$\\frac{ep_{a}}{ep_{b}}=%.1fep_{a}$' % (1/epb))

        plt.annotate('$%.1f$' % (epb), xy=(elist[-1]-0.05, KallList[-1]-0.03), color=colorlist[j % len(colorlist)], size=9)
        j += 1
        
        plt.savefig("./pic/"+'ep%d_%d.png'%(epa*10,epb*10), dpi=400, bbox_inches='tight')

    

    plt.xlabel('${\\xi}$')
    plt.ylabel('${K^{LS}/K_{\max}}$/(N/(mm${\\cdot\\mu}$m))')
    #plt.legend()

    #plt.savefig("./pic/"+'epdiv%d.png'%(epa*10), dpi=400, bbox_inches='tight')
    plt.show()
    '''
    
